package scu.maqiang.possion;

import scu.maqiang.fes.*;
import scu.maqiang.mesh.*;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.ParamCheck;
import scu.maqiang.numeric.SRMatrix;

public class Possion2DQ41 {

	public Possion2DQ41(Mesh2Q4 mesh) {
		this.mesh = mesh;
		fs = new FES2Q41(mesh);
		A = new SRMatrix(fs.GetNdof());
		x = new double[fs.GetNdof()];
		gradient = new double[2][fs.GetNdof()];
		RHS = new double[fs.GetNdof()];
	}
	
	public static void main(String[] args) {
        int N = 8, label = 10;
        double ep = 1.0 / N;
       
//        //使用lambda函数来指标材料系数分布
//        RegionFunc regionFunc = (xy, param) -> {
//            double xCell = xy[0] * N - (int)(xy[0] * N);
//            double yCell = xy[1] * N - (int)(xy[1] * N);
//            return (Math.abs(xCell - 0.5) < 0.25 && Math.abs(yCell - 0.5) < 0.25)? 1.0:-1.0;
//        };
//        
//		Mesh2Q4Generator mg = new Mesh2Q4Generator();
//		Mesh2DQ4 mesh2D = mg.square2D(100, 100);
//        mesh2D.setDomainLabel(regionFunc, null, label);
//        
//        FES2Q41 fs = new FES2Q41(mesh2D);
//        SRMatrix A = new SRMatrix(fs.GetNdof());
//
//        GeneralCoefFunc regionCoef = (xyz, relabel, param) -> relabel == label ? 0.01 : 1.0;
//        GeneralCoefFunc ConstCoef = (xyz, relabel, param) -> 1.0e0;
//        GeneralCoefFunc[] gcf = {regionCoef};
//        gcf[0] = ConstCoef;
//        
//        GeneralCoefFunc coefFunc = (xyz, relabel, param) -> 2 * xyz[0] * (1 - xyz[0]) + 2 * xyz[1] * (1 - xyz[1]);
//        fs.assembleHeatStiff(gcf, null, BVPType.COMMON, A);
//        //fs.assembleHeatStiff(new double[] {1.0}, BVPType.COMMON, A);
////        System.out.println(A);
////        System.out.println(A.makeEntrySymmetric());
//        double[] RHS = new double[fs.GetNdof()];
//        //GeneralCoefFunc ConstCoef = (xyz, relabel, param) -> 1.0;
////        coefFunc = (xy, relabel, param) ->  {
////        	double a = 400.0;
////        	double r = Math.sqrt((xy[0] - 0.5) * (xy[0] - 0.5) + (xy[1] - 0.50) * (xy[1] - 0.5));
////        	return 4 * 400 * 400 * (1 - a * r * r) * Math.exp(-400.0 * r * r);
////        };
//        gcf[0] = coefFunc;
//        fs.assembleHeatSource(gcf, null, BVPType.COMMON, RHS);
//        A.makeEntrySymmetric();
//        fs.applyBCZero(A, RHS, Direct.X, BCType.RCE, 1, 2, 3, 4);
//        
////        System.out.println(A);
//
//        double[] x = new double[fs.GetNdof()];
//        
//        //IterSSolver solver = new IterSSolver(1.0e-15, 1000000);
//        
//        //CSRMatrix mat = new CSRMatrix(A);
//        //solver.PCGSSOR(A, RHS, x, 1.5e0, 1);
//        //Arrays.fill(x, 0.0);
//        //solver.PCGSSOR_2(A, RHS, x, 1.5e0, 1);
//        
//        //System.out.println(A);
//        
//        //System.out.println(mat);
//        
//        //Arrays.fill(x, 0.0);
//        //solver.CG(mat, RHS, x, 1);
//        
//        //Arrays.fill(x, 0.0);
//        //solver.CG(A, RHS, x, 1);
//        
//        Arrays.fill(x, 0.0);
//        NewIterSSolver newSolver = new NewIterSSolver(A);
//        newSolver.PCGSSOR(RHS, x, 1.5e0, 1);
//
//        double[][] gradX = new double[2][fs.GetNdof()];
//        fs.computeGradient2(x, gradX);
//        mesh2D.toTecplot("Solutions.dat", x);
//        mesh2D.toTecplot("SolutionsGradient.dat", gradX);
//        
//        System.out.println("Gradient: ");
//        System.out.println("Dx min: " + MVO.min(gradX[0]) + "\tmax: " + MVO.max(gradX[0]));
//        System.out.println("Dy min: " + MVO.min(gradX[1]) + "\tmax: " + MVO.max(gradX[1]));

		Mesh2Q4 mesh2D = new Mesh2Q4().square2D(50, 50);

		ScalarFunc ConstFunc = (xyz, relabel, param) -> 1.0e0;
		ScalarFunc SourceFunc = (xyz, relabel, param) -> 2 * xyz[0] * (1 - xyz[0]) + 2 * xyz[1] * (1 - xyz[1]);
		Possion2DQ41 p2DT31 = new Possion2DQ41(mesh2D);
		p2DT31.DirichletProblem(BVPType.COMMON, ConstFunc, null, SourceFunc, null, null, null, null, "MyResultQ4.dat");
	}
	
	/**
	 * Possion方程Dirichlet边界条件求解器
	 * <p>
	 * &part;<sub>i</sub>(a<sub>ij</sub>&part;<sub>j</sub>u) = f &#915;
	 * 
	 * @param kxxFunc
	 * @param paramKxx
	 * @param sourceFunc
	 * @param paramSource
	 * @param tp
	 * @param DirichletBoundaryLabel
	 * @param fixedFunc
	 * @param paramFixed
	 * @param reportFile
	 */
	public void DirichletProblem(BVPType tp, ScalarFunc kxxFunc, double[] paramKxx, ScalarFunc sourceFunc,
			double[] paramSource, int[] DirichletBoundaryLabel, ScalarFunc fixedFunc, double[] paramFixed,
			String reportFile) {
		fs.assembleStiff(new ScalarFunc[] { kxxFunc }, paramKxx, tp, A);
		if (sourceFunc != null) {
			fs.assembleSource(new ScalarFunc[] { sourceFunc }, paramSource, tp, RHS);
		}
		fs.applyBC_MBN(A, Direct.All, DirichletBoundaryLabel);

		if (fixedFunc == null) {
			fs.applyBC_MBN(RHS, Direct.All, 0.0, DirichletBoundaryLabel);
		} else {
			fs.applyBC_MBN(RHS, Direct.All, new ScalarFunc[] { fixedFunc }, paramFixed, DirichletBoundaryLabel);
		}
		solveAndComputeGradientAndExport(reportFile);
	}

	public void DirichletProblem(BVPType tp, double kxx, double source, int[] DirichletBoundaryLabel,
			double fixedBV, String reportFile) {
		fs.assembleStiff(new double[] { kxx }, tp, A);
		fs.assembleSource(new double[] { source }, tp, RHS);
		fs.applyBC_MBN(A, Direct.All, DirichletBoundaryLabel);
		fs.applyBC_MBN(RHS, Direct.X, fixedBV, DirichletBoundaryLabel);
		solveAndComputeGradientAndExport(reportFile);
	}

	private void solveAndComputeGradientAndExport(String reportFile) {
		NewIterSSolver solver = new NewIterSSolver(A, true);
		solver.PCGSSOR(RHS, x, 1.50, 1);
		//DirectSSolver solver = new DirectSSolver(A, fs.getnBand());
		//solver.PLUSolve(RHS, x, 1);
		fs.computeGradient2(x, gradient);
		System.out.println("Solution: ");
		System.out.println("x0 min:" + MVO.min(x) + "\t" + "max: " + MVO.max(x));
		System.out.println("Solution Gradient: ");
		System.out.println("Dx min:" + MVO.min(gradient[0]) + "\t" + "max: " + MVO.max(gradient[0]));
		System.out.println("Dy min:" + MVO.min(gradient[1]) + "\t" + "max: " + MVO.max(gradient[1]));
		double[][] mat = new double[][] { x, gradient[0], gradient[1] };
		mesh.toTecplot(reportFile, mat);
	}

	public void GeneralPossionProblem(BVPType tp, ScalarFunc kxxFunc, double[] paramKxx, ScalarFunc sourceFunc,
			double[] paramSource, int[] DirichletBoundaryLabel, ScalarFunc fixedFunc, double[] paramFixed,
			int[] NeumannBoundaryLabel, ScalarFunc fluxFunc, double[] paramFlux, int[] RobinBoundaryLabel,
			ScalarFunc robinFluxFunc, double[] paramRobinFlux, ScalarFunc convectFunc,
			double[] paramConvectFunc, String reportFile) {

		fs.assembleStiff(new ScalarFunc[] { kxxFunc }, paramKxx, tp, A);
		if (sourceFunc != null) {
			fs.assembleSource(new ScalarFunc[] { sourceFunc }, paramSource, tp, RHS);
		}
		if (NeumannBoundaryLabel != null) {
			fs.assembleFlux(new ScalarFunc[] { fluxFunc }, paramFlux, tp, RHS, NeumannBoundaryLabel);
		}
		if (RobinBoundaryLabel != null) {
			fs.assembleBoundaryMass(new ScalarFunc[] { convectFunc }, paramConvectFunc, tp, A,
					RobinBoundaryLabel);
			fs.assembleFlux(new ScalarFunc[] { robinFluxFunc }, paramRobinFlux, tp, RHS, RobinBoundaryLabel);
		}

		fs.applyBC_MBN(A, Direct.All, DirichletBoundaryLabel);
		fs.applyBC_MBN(RHS, Direct.X, new ScalarFunc[] { fixedFunc }, paramFixed, DirichletBoundaryLabel);
		solveAndComputeGradientAndExport(reportFile);

	}

	public void GeneralPossionProblem(BVPType tp, double kxx, double source, int[] DirichletBoundaryLabel,
			double fixedBV, int[] NeumannBoundaryLabel, double fluxValue, int[] RobinBoundaryLabel,
			double robinFluxValue, double convectValue, String reportFile) {
		fs.assembleStiff(new double[] { kxx }, tp, A);
		fs.assembleSource(new double[] { source }, tp, RHS);
		if (NeumannBoundaryLabel != null) {
			fs.assembleFlux(new double[] { fluxValue }, tp, RHS, NeumannBoundaryLabel);
		}
		if (RobinBoundaryLabel != null) {
			fs.assembleBoundaryMass(new double[] { convectValue }, tp, A, RobinBoundaryLabel);
			fs.assembleFlux(new double[] { robinFluxValue }, tp, RHS, RobinBoundaryLabel);
		}
		
		A.sort();
		System.out.println(A);
		fs.applyBC_MBN(A, Direct.All, DirichletBoundaryLabel);
		fs.applyBC_MBN(RHS, Direct.X, fixedBV, DirichletBoundaryLabel);
		solveAndComputeGradientAndExport(reportFile);
	}
	
	public void GeneralPossionProblem(BVPType tp, double kxx, double source, int[][] DirichletBoundaryLabel,
			double[] fixedBV, int[][] NeumannBoundaryLabel, double[] fluxValue, int[][] RobinBoundaryLabel,
			double[] robinFluxValue, double[] convectValue, String reportFile) {
		fs.assembleStiff(new double[] { kxx }, tp, A);
		fs.assembleSource(new double[] { source }, tp, RHS);
		double[] value = {0.0};
		if (NeumannBoundaryLabel != null) {
			int n = fluxValue.length;
			ParamCheck.checkEqual(n, NeumannBoundaryLabel.length);
			for(int i = 0; i < n; i++) {
				value[0] = fluxValue[i];
				fs.assembleFlux(value, tp, RHS, NeumannBoundaryLabel[i]);
			}
			
		}
		if (RobinBoundaryLabel != null) {
			int n = robinFluxValue.length;
			ParamCheck.checkEqual(n, convectValue.length);
			ParamCheck.checkEqual(n, RobinBoundaryLabel.length);
			for(int i = 0; i < n; i++) {
				value[0] = convectValue[i];
				fs.assembleBoundaryMass(value, tp, A, RobinBoundaryLabel[i]);
				value[0] = robinFluxValue[i];
				fs.assembleFlux(value, tp, RHS, RobinBoundaryLabel[i]);
			}

		}
		
		int n = DirichletBoundaryLabel.length;
		ParamCheck.checkEqual(n, fixedBV.length);
		for(int i = 0; i < n; i++) {
			fs.applyBC_MBN(A, Direct.All, DirichletBoundaryLabel[i]);
			fs.applyBC_MBN(RHS, Direct.X, fixedBV[i], DirichletBoundaryLabel[i]);
		}
		
		A.sort();
		System.out.println(A);
		solveAndComputeGradientAndExport(reportFile);
	}
	
	public double[] exportSolution() {
		return x.clone();
	}
	
	public double[][] exportGradient() {
		double[][] gd = new double[2][];
		gd[0] = gradient[0].clone();
		gd[1] = gradient[1].clone();
		return gd;
	}
	
	Mesh2Q4 mesh;
	FES2Q41 fs;
	SRMatrix A;
	double[] RHS;
	double[] x;
	double[][] gradient;


}
